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Abstract. I present an algorithm for inverting the luminosity function for white 
dwarfs to obtain a maximum likelihood estimate of the star formation rate of the host 
stellar population. The algorithm is of the general class of Expectation Maximization, 
and involves iteratively improving an initial guess of the star formation rate. Tests show 
that the inversion results are quite sensitive to the assumed metallicity and initial mass 
function, but relatively insensitive to the initial-final mass relation and ratio of H/He 
atmosphere white dwarfs. Application to two independent determinations of the Solar 
i-C ■ neighbourhood white dwarf luminosity function gives similar results: the star formation 

rate is characterised by an early burst, and more recent peak at 2-3 Gyr in the past. 



1. Introduction 

The white dwarf luminosity function (WDLF) is a useful tool for determining the age 
of a population of stars. The magnitude at which the function terminates is highly 
time-dependent, and by fitting the faint end with theoretical WDLF models of different 
ages one can obtain a statistical estimate of the age of the population without having 
to determine the total age of any individual white dwarf, which is considerably more 
difficult. This technique has been applied both to sin gle burst populations such as open 
clusters (IBedin et al. l2010l:lG arcfa-Berro et al.ll20 10) and continuous populations such 
as the Galactic disk (IQswalt et al.lll996l ; lKnox et al.lll999h . 

T he standard equation for modelling the WDL F for a given star formation history 



is (e.g. llben & Laughlinll 19891 : Fontaine et al.H2001h 



3>(M bol )= f J^^To-tcooi-tusJtMdM (1) 
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where ®(Mt, i) is the number density of WDs at magnitude M bo] . The derivative inside 
the integral is the characteristic cooling time for WDs, is the star formation rate 
(SFR) at time t and <p is the initial mass function (IMF). The integral also depends on 
the lifetimes of main sequence progenitors as a function of mass and metallicity ?ms> 
the WD cooling times as a function of mass and luminosity f coo i, the initial-final mass 
relation m(M) and the total time since the onset of star formation Tq. The integral is 
over all main sequence masses that have had time to produce WDs at the present day, 
with the magnitude-dependent lower limit corresponding to the solution of 
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and the upper limit M u « 7 M Q . 

From various studies of this equation ( Iben & Laughlin 1989; iNoh & Scalolll990h 
it is known that the faint end of the WDLF is mostly insensitive to the SFR, and is 
determined mainly by the total age of the population To. WDs at these magnitudes 
are uniformly old, and are the remains of high mass main sequence stars that formed 
right at the onset of star formation. It is for this reason that the faint end provides 
the most constraint on the total age. The picture is considerably more complicated at 
brighter magnitudes, because the WDs are a mixture of ages: both young, high mass 
WDs that are produced by recently-formed MS progenitors, and old, low mass WDs 
that are produced b y low mass MS stars that formed at early times. It was found by 
Noh & Scald (Il990h that time variations in the SFR may leave imprints in the WDLF 



at these magnitudes, and by forward modelling methods they interpreted a marginal 
feature in the WDL F at Mh n i ~ 1 as e vidence for a burst of star formation 0.3 Gyr 
ago. According to INoh & Scald (1990) and equation \T\ the shape of the WDLF at 
intermediate magnitudes is also strongly affected by the cooling rates of WDs, and 
it is possible that features in the WDLF may be i nterpreted as evidence of additional 
WD cooling mechanisms (see e.g. Ilsern et alj|2008l and the contribution of Miguel & 
Bertolami to these proceedings). 

This paper presents results of ongoing work on a strategy to invert the WDLF 
to obtain a direct estimate of the time varying SFR. This work is driven by two related 
questions: given current WD cooling models, what constraint can features in the WDLF 
(at all magnitudes) place on the time varying SFR? And as a corrolary to this: can 
features in the WDLF be explained exclusively by time variations in the SFR, or are 
additional cooling mechanisms required? 



2. White Dwarf Luminosity Function Inversion Algorithm 

To a first approximation, the two parameters that determine the total age of a WD are 
the present day bolometric magnitude, and the mass. These can be used to determine 
both the total WD cooling time and the time spent on the main sequence. The approach 
to inverting the WDLF presented here is based on the observation that if the distribution 
of WD mass was known at all magnitudes, then the WDLF could be immediately trans- 
formed to the SFR. As this quantity is generally not known observationally, this direct 
approach can't be used. Instead, we use the inversion techniq ue known as Expectation 
Maximization (Demps ter et al.l [l977l : iDo & B atzogloul l2008h . which is used to obtain 



maximum likelihood estimates of the solution to inverse problems in the presence of 
missing data. 

This approach involves iteratively refining an initial guess of the SFR. The general 
procedure for each iteration is as follows. The starting point is an initial guess of the 
star formation rate tfro, 

<Ao = Mt), (2) 

where in the present work \]/q is flat, i.e. a constant star formation rate. This is combined 
with the initial mass function (p to get the joint mass and formation time distribution of 
main sequence progenitors Pms, where 



PM$(M MS ,t) = (/>(Mms)iAo(0 



(3) 
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Using standard rules of probability density functions, we can transform this to the joint 
mass and bolometric magnitude distribution of white dwarfs Pwd at the present day: 



Pwd(M W d, Afboi) - Pms(M M s, 0- 



d(M MS , t) 



d(M WD , Afboi) 



(4) 



This function can be seperated into a product of the marginal luminosity distribution 
and the mass distribution conditioned on luminosity, 

P W d(M wd , Mboi) - O sim (M bol )PwD(M W D|Mboi) (5) 

The quantity O s j m is just the WDLF for the initial guess SFR model, up to a normalisa- 
tion factor. The next crucial step is to replace this with the observed WDLF ® t, s to get 
the updated WD distribution P' WD - 

P' WD (M WD ,M bo i) = b s (Mboi)PwD(M WD |M b oi) (6) 

This updated WD distribution has the same marginal luminosity distribution as the ob- 
served WDLF, and the magnitude-dependent mass distribution derived from the initial 
guess star formation rate model. We can now invert this distribution to obtain the up- 
dated distribution for main sequence stars P' MS again using standard transformation 
rules: 

<9(M WD , Mboi) 

P MS ( m ms, t) = P WD (M WD , Mboi) 



(7) 

d(M MS ,t) 

The final step is to marginalise P' MS over the main sequence mass, to obtain the updated 
star formation rate model if/\ : 



1 r M ™ 

1 - A{t) jM .fetime/,) 



(8) 



The integral is over all main sequence stars that produce WDs at the present day. In 
the present work, the upper limit - 7M Q , and the variable lower limit M l ^f ime (t) 
correponds to the mass of the main sequence star with lifetime t. The factor A corrects 
for low mass MS stars that don't form WDs at the present day, and is calculated 

A(f)= 4>(Mus)dM MS (9) 

Jm ms 

where the lower mass limit in this case is set to M?S? - 0.6 M Q . The star formation rate 
recovered by this algorithm therefore only represents stars more massive than 0.6 M Q . 
It is also non-parametric, in the sense that it does not enforce any particular functional 
form on the recovered SFR. 



2.1. White Dwarf Atmosphere Types 

Along with the mass and present luminosity, the H/He atmosphere type is a third param- 
eter affecting the total age of a WD. This has a significant effect at larger cooling ages 
(> 6 Gyr depending on the choice of models), with H atmosphere WDs being brighter 
at a given cooling age. The two atmosphere types can be included in the algorithm in a 
relatively straightforward manner. We calculate Pwd in equation |4] separately for each 
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atmosphere type, then take a linear combination to obtain the total P WD for the mixed 
atmosphere population, where 

Py fD = aP^ D + (l-a)P^ 3 (10) 

The factor a fixes the relative abundance of H and He WDs at birth, though their ratio 
changes with luminosity due to the two types cooling at different rates. A value of 
a = 0.5 is used in the present work. 



3. Validation with Synthetic Data 



In order to test the accuracy of the recovered SFR, we have generated a set of synthetic 
WDLFs using a range of different known input SFR models. This allows us to check the 
performance of the algorithm in tightly controlled noise conditions, and the sensitivity 
to uncertainties in the various modelling inp uts. I n this work we use the WD cooling se- 
quences described in iTremblay et all (1201 lh and iBergeron et all (1201 lh and references 
therein (see also www.astro.umontreal.ca/~bergeron/CoolingModels). Figure 
Q] shows the results from two tests in noise-free conditions, designed as a proof of con- 
cept to verify that the algorithm works in principle on realistic models of the SFR. In 
each case, the synthetic WDLF (not shown) has a magnitude binning of AMb i = 0.5, 
chosen to match the observed WDLF resolution in recent studies. The algorithm per- 
forms well on smoothly varying SFRs like the exponential decay model (left). The 
overall form of the fractal SFR model on the right is recovered by the inversion, al- 
though at older times high frequency components in the underlying SFR are lost and 
the algorithm only measures a moving average. This is a fundamental limit of the algo- 
rithm arising from the finite magnitude resolution in the WDLF. 




Lookback time [Gyr] Lookback time [Gyr] 

Figure 1 . Testing convergence of inversion algorithm using synthetic WDLF data. 
The black lines show the underlying SFR model used to generate the artificial WDLF 
in each case; the red lines show the recovered SFR model on the final iteration of the 
algorithm. In these tests, the algorithm converged in 27 and 32 steps. 



Inverse problems are notoriously sensitive to noise. We have carried out a com- 
prehensive campaign of tests to assess the effect of observational errors and modelling 
uncertainties on the algorithm performance. To summarise, observational errors on the 
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scale of recent WDLF measurements in the Solar neighbourhood do not catastrophi- 
cally degrade the performance of the algorithm, and we are able to correctly estimate 
the true error on the inverted SFR. In terms of modelling parameters, the algorithm is 
most sensitive to uncertainties in the IMF and progenitor metallicity, and less sensitive 
to the initial-final mass relation and H/He atmosphere ratio. The full results will be 
published separately. 



4. The Solar Neighbourhood 



This algorithm has been applied to two recent meas urem ents of the WDLF fo r the Solar 
neighbourhood: that of Rowell & Hambly <l201lh and lHarris et all ((2006) (hereafter 
RH11 and H06). The results are shown in figure |2j Both SFRs show a similar form, 
being characterised by an early burst and a more recent peak at 2-3 Gyr in the past. 
The difference in magnitude is due to the significant incompleteness (~50%) of the 
RH11 sample with respect to that of H06. In both cases, the maximum lookback time 
is fixed at 9 Gyr. In figure [3] the best fit synthetic WDLFs found on convergence of 
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Figure 2. Star formation rates recovered from two recent determinations of the 
Solar neighbourhood WDLF, that of RH1 1 and H06. The filled regions show the l<x 
uncertainty. In these tests, the algorithm converged in 1 1 and 28 steps. 



the algorithm are plotted over the observed WDLFs. The inset panels show the ratio of 
the two functions. The WDLFs are fitted very well by the algorithm, particularly in the 
case of the H06 WDLF, and there appears to be no significant over- or under-abundance 
of WDs that remains unaccounted for. 



5. Conclusion 

We have presented preliminary results of work to invert the WDLF. This represents a 
new method of analysing the star formation history of the Solar neighbourhood and 
WD populations more generally, one that is essentially independent of existin g in- 
version methods tha t use main sequence stars, such as iHernandez et all (l2000l) and 
ICignoni et all ( 20061) . Application of the algorithm to the Solar neighbourhood WDLF 
yields a SFR characterised by an early burst and a recent (~ 2 - 3 Gyr ago) peak. Future 



6 



N. Rowell 




Figure 3. Best fit WDLFs obtained from converged SFR models show a good fit 
to the observed WDLF in each case. 



development work will include making the maximum lookback time a free parameter. 
We also plan to compare results for different sets of WD cooling models, which may 
turn out to be the largest uncertainty in the recovered SFR. It would also be interest- 
ing to apply the method to other WD populations such as the thick disk, spheroid and 
clusters. Single burst populations in particular would provide a useful benchmark. 
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